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ABSTRACT 


The design of the nozzle-afterbody section for a hypersonic transport such as the NASP is 
currently underway, and is being conducted using both computational fluid dynamics and cold, 
non-reacting, simulant gas experimental models. In this study, computations are performed for 
a cold gas simulation of a scramjet afterbody flowfield and compared with the results obtained 
from an experimental study of scramjet module nozzle-afterbody flows. The expansion of 
a supersonic flow through an intemal/extemal nozzle-afterbody configuration and its viscous 
mixing with a hypersonic freestream flow of air is computed using two and three-dimensional 
upwind, finite volume, Navier-Stokes schemes. The Reynolds stresses are represented by a 
Baldwin-Lomax algebraic turbulence model with modifications to account for separated flow, 
multiple wall geometry, and turbulent wake flow. Two-dimensional computations for over- 
expanded (off design) and under-expanded nozzle flows, are performed on flow adapted grids, 
and three-dimensional results are obtained for a half-span nozzle model on a fixed grid. The 
results obtained from the adapted grid computations show improved accuracy and resolution. The 
computed pressure distributions compare favorably with experimental results. Furthermore, the 
results demonstrate that the solutions obtained from the computational fluid dynamics algorithms 
used in this study, can be used to expand the database for these types of nozzle-afterbody 
configurations. 
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LIST OF SYMBOLS 


a = local speed of sound 
C p = coefficient of pressure 
C v = constant volume specific heat 
e t = total energy per unit volume 
e = internal energy per unit volume 
F, G, H = inviscid flux vectors 
F y , G v , H v = viscous flux vectors 
J = transformation Jacobian 

k = coefficient of thermal conductivity; also tension spring constant 
for adaptive grid algorithm 
M = Mach number 
n = normal distance from a boundary 
p = pressure 
Pr = Prandd number 
q = heat flux 

Q = vector of conserved variables representing mass, momentum, 
and total energy per unit volume 
R = universal gas constant 
Re = Reynolds number 
t = time 

T = static temperature 
u,v,w = cartesian velocities 

U, V, W = contravariant velocities; V also denotes the volume of a computational cell 
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GREEK 


6 = boundary layer thickness; also central difference operator 

7 = ratio of specific heats 

r), £ = curvilinear coordinates 
A = bulk viscosity 
r = shear stress 
p = density 
<x>= vorticity 

p = effective viscosity (pi + p x ) 
pi = molecular viscosity 
Pi = eddy (turbulent) viscosity 

SUBSCRIPTS AND SUPERSCRIPTS 

i, j, k = spatial indices 
n = time level 
w = wall condition 
oo = freestream conditions 
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Chapter 1 
INTRODUCTION 


General background information and motivation for this study are presented in this chapter. 
A brief survey of the related current literature, both experimental and computational, is provided 
along with the specific objectives of this study 

1.1 Background 

The recent resurgence of interest in hypersonic aerodynamics has come about largely in 
part due to the development of the National Aerospace Plane (NASP). The design of this type 
of aircraft will rely heavily on the use of computational fluid dynamics, since the operating 
conditions (Mach number, Reynolds number, enthalpy levels, etc.) prohibit the use of most of 
the conventional experimental facilities to obtain the required data for design analysis. 

One of the major design tasks involved in the development of a hypersonic air-breathing 
aircraft is the integration of the engine and the airframe. This is necessary in order to reduce 
excessive drag and weight due to the high Mach numbers at which the aircraft will be traveling. 
Preliminary designs for an integrated propulsion system include a forebody/inlet system in which 
the forebody of the aircraft is used to partially compress the air before it enters the engine module. 
This is shown in the sketch of Fig. 1.1. The oncoming air is further compressed by a series of 
wedges and struts, and then mixed with a fuel and ignited in the combustor region. Finally, the 
high pressure combustion products are expanded through the combustor exit nozzle and over the 
airframe afterbody configuration. The overall propulsive efficiency of the nozzle is determined, 
to a large extent, by the exhaust plume flow over this afterbody section. This study is concerned 
primarily with this nozzle-afterbody region. 

The design and testing of a scramjet nozzle-afterbody section using actual engine combustion 
products is impractical in a conventional wind tunnel. The actual chemistry and high total 
enthalpy levels of the exhaust products would be quite difficult to match in a scaled test 


section. However, several alternatives do exist. A simulant gas can be substituted for the 
actual combustion products, provided that dynamic and thermodynamic similitude are enforced. 
Perhaps a more economical alternative would be to do the preliminary design analysis using 
computational fluid dynamics (CFD). Since there is currently very little experimental data for 
very high Mach number flows, some means of calibrating and validating these CFD codes must 
be achieved before they can be used with complete confidence in this design process. 

1,2 Computational and Experimental Literature on Nozzle-Afterbody Flows 

In the 1970’s, a study was undertaken by Grumman Aerospace Corporation under a NASA 
contract to develop an experimental cold gas simulation technique for scramjet exhaust flows 
[1-3]*. The purpose of this study was to define a method for accurately simulating pressure 
distributions on the nozzle-afterbody surfaces of a hypersonic scramjet aircraft using a cold 
substitute simulant gas, as opposed to the hot combustion product gases from the scramjet 
engine. It was determined in this study, that in addition to the usual nondimensional similitude 
parameter requirements for inviscid flows (i.e. Mach numbers, pressure ratios, temperature 
ratios, etc.), that the ratio of specific heats (7) of the combustion products must also be matched 
by the simulant gases. This concept is discussed in further detail in Appendix A. It was also 
determined in this study that the surface pressures were relatively insensitive to small changes 
in the thermodynamic properties of the gases, but are very sensitive to flow perturbations caused 
by the nozzle geometry. 

An extension of this work was carried out at the NASA Langley Research Center by Cubbage 
et al. [4—5] and Pittman [6]. Experimental data was obtained for a scaled scramjet nozzle- 
afterbody flowfield using both air and a Freon/Argon mixture as the simulant gas. Static pressures 
were measured on the afterbody surface, for both two-dimensional and three-dimensional flows, 
with various nozzle-afterbody geometries. The experimental data of [5] was compared with 
numerical data from two-dimensional Navier-Stokes and Euler algorithms. It was concluded 

* The numbers in brackets indicate references 
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from the experimental results of [6] that the pressure drop in a supersonic expansion will be 
smaller for the Freon/Argon simulant gas than for air, and hence results in greater overall 
afterbody forces. 

Novak and Cornelius [7] conducted subsonic and transonic experiments on a similar af- 
terbody model using air as the simulant gas. Their work included three-dimensional laser ve- 
locimetry measurements for the subsonic test portion of the nozzle-afterbody flow field. Very 
little experimental data exists for off-surface flowfields of supersonic/hypersonic mixing flows. 
This data would be very important in calculating exhaust plumes and shear layers, which may 
cause interference on aerodynamic control surfaces of the aircraft. 

Extensive numerical studies have been conducted on rocket and nozzle based flows without 
solid afterbody surface interaction, that is, the flow from a nozzle exit plane which flows directly 
into the freestream. Most of these studies focus mainly on the thrust and thrust vectors generated 
by the nozzle itself, and not on the actual afterbody flow field plume structures. However, there 
have been several of these types of studies which did attempt to analyze the downstream flow. 

Deiwert [8] used the thin- shear-layer formulation of the Navier-Stokes equations in the study 
of supersonic axisymmetric flow over boattails containing a centered propulsive jet. Solutions 
were presented for jet flows expanding supersonically into a low pressure supersonic freestream 
flow, with an in-depth analysis of the afterbody flow structure. Comparison with experimental 
data showed good agreement in many of the key flow features such as exhaust plume shape 
and structure, and the location of the external compression shocks. However, the quantitative 
comparison of nozzle flow exit angle was poor, and thought to be due to the improper modeling 
of turbulent transport phenomena in the region of separated flow at the base of the nozzle. 

Hoffman et al. [9] computed the afterbody flowfield of an axisymmetric rocket nozzle 
flowfield using a complete Navier-Stokes formulation algorithm. They presented the results of 
a grid refinement study which showed that the grid resolution is very important in regions of 
shear layers and recirculation. The computed results obtained on a “quality” grid, where the grid 
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point spacing was very dense in regions where large gradients were expected to occur, compared 
much more favorably with the experimental results. 

Goldberg et al. [10] studied afterbody flowfields at sonic and supersonic Mach numbers. 
They presented results obtained from an upwind Navier-Stokes solver using alternatively a k-e 
or a Baldwin-Lomax turbulence model. Comparisons with experimentally obtained pressure data 
were very good, however both turbulence models were deemed inadequate near large regions of 
flow separation. They also found that finer grid clustering and more grid points were needed in 
order to resolve the intricate details of the flow 

In recent years, with the renewed interest in hypersonics and propulsion airframe integration, 
there have been numerous computational studies related to nozzle-afterbody flows. Barber and 
Cox [11] presented an overview of computational works being conducted towards the design 
of hypersonic airbreathing aircraft, with a section specifically devoted to the integration of the 
propulsion system. Povinelli [12] provided a summary of the current computational works which 
are directly related to the propulsion systems, including the nozzle-afterbody section. Many of 
the works have concentrated on one particular aspect of the flow or its approximation, such as 
the turbulence modeling, adaptive gridding schemes, shear layer analysis, or real gas effects. 
Several of the studies which are directly related to this work are given here. 

Ray et al. [13] presented two and three-dimensional results based on Euler (inviscid) 
calculations of single and multiple module scramjet afterbody flows using a Freon/Argon simulant 
gas mixture. Their results show good agreement with the experimental data of [3], however their 
calculations are based on the nozzle flow expansion into quiescent air, not the actual hypersonic 
flow of air. Also, since they are treating the flow as inviscid, shear layers and regions of 
separated flow can not be computed. 

Hsu [14] conducted a study of adaptive gridding refinements based on various nozzle- 
afterbody configurations. Several cases were calculated with and without the use of adaptive 
grids, and the end result in most cases, especially those which contained free shear layers and 
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strong shocks within the flow, was the fact that the use of flow adaptive grids is essential in 
improving solution accuracy. 

Baysal et al. [ 15 ] performed a two-dimensional analysis of the experimental work carried 
out by Cubbage et al. [ 4 - 5 ] and Pittman [6]. In the work of [ 15 ] the flow was analyzed using 
two different algorithms. The first code used was an implicit, upwind, finite volume Navier- 
Stokes solver which assumed a constant specific-heat-ratio (7). The second code was an explicit 
MacCormack based Navier-Stokes solver, which included coupled species continuity equations 
to account for variable 7 due to the mixing of the simulant gas (Freon/Argon mixture) with the 
external freestream air. This work was extended by Baysal et al. [ 16 , 17 ] to include three- 
dimensional flow and adaptive gridding methodology in the analysis. The two-dimensional 
computational meshes were refined using flow adaptive grids to enhance the computational 
solution in the regions of large gradients, and to reduce the overall computational error. Three- 
dimensional flow solutions were computed for the nozzle-afterbody test section and compared 
with the three-dimensional experimental surface pressure data of [ 5 ]. 

Harloff et al. [ 18 ] conducted two-dimensional Navier-Stokes analysis of scramjet nozzle 
flow fields at design and off-design conditions. In this study, both the nozzle exhaust flow 
and the external flow were assumed to be air. Nozzle afterbody flowfields were computed over 
a range of nozzle exit Mach numbers for over-expanded and under-expanded flows. Nozzle 
efficiencies were computed for all of the cases, however no experimental data was given for 
comparison with numerical results. 

Edwards [ 19 ] studied the exhaust plume/afterbody interactions using the thin-layer Navier- 
Stokes assumption with a coupled species continuity equation. This allowed for the solution of 
a binary gas mixture flow. Computations were performed for the external flow of air (7=1.4) 
mixing with a simulant gas (7=1.26). Additional computations were done with the simulant gas 
assumed to be air (7=1.4). It was concluded, from the computational results and a kinetic theory 
of gases rationalization, that the pressure drop in supersonic expansion will be smaller for the 
7=1.26 gas than for the 7=1.4 gas, and hence produced greater afterbody forces. 
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Ruffin et al. [20] have used two-dimensional and three-dimensional upwind Navier-Stokes 
solvers to do a preliminary analysis of a planned nozzle-afterbody experiment. They conducted 
two-dimensional parametric studies over various Mach numbers, pressure ratios, and afterbody 
ramp angles to help determine the experimental model loads and optimum afterbody ramp angle 
and length. Three-dimensional calculations were performed to predict the shape of the jet plume 
and the flow spillage from the windward side of the model into the expanding flow region. The 
three-dimensional results were also to be used to determine optimal locations for experimental 
probes and flow measurement devices, and to aid in the design of side flow fences which are 
used to minimize the flow spillage. 

1.3 Objectives 

The main emphasis of this research effort is focused on developing a further understanding of 
scramjet nozzle-afterbody flowfields through the use of computational fluid dynamics. This study 
is being conducted in parallel with wind tunnel tests of scramjet nozzle-afterbody flowfields, so 
that the computational results obtained can be compared with experimental results in order to 
benchmark the solution algorithms. 

The specific objectives of this study include: 

(a) developing a means of computationally analyzing the two-dimensional viscous mixing of a 
scramjet exhaust flow with that of an external freestream flow. 

(b) applying an adaptive gridding methodology to the results obtained in (a) to enhance the 
solution and reduce the computational error. 

(c) extending the two-dimensional analysis to three-dimensions, with the intention of using the 
computational methods to augment the data base on scramjet nozzle-afterbody flows. 


The basic physical characteristics of the scramjet nozzle-afterboby flowfields are discussed 
in Chap. 2. The governing equations and corresponding boundary conditions are given in Chap. 
3, and the basic algorithm used in the solution of these equations, along with the grid generation 
methods are detailed in Chap. 4. Results and discussion are presented in Chap. 5, and some 
conclusions and recommendations are given in Chap. 6. 
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Chapter 2 

PHYSICAL DESCRIPTION 

Propulsion airframe integration is taking a vital role in the development of a hypersonic 
airbreathing aircraft such as the National Aerospace Plane. One of the underlying concepts 
of propulsion airframe integration is the incorporation of the external nozzle modules with the 
airframe controlling surfaces. In this study, the afterbody surface of the fuselage acts to expand 
the supersonic turbulent flow scramjet exhaust, hence becoming part of the nozzle. The expansion 
of exhaust products over the external nozzle-afterbody surface will produce large forces and 
pitching moments on the aircraft, therefore some means must be developed to accurately predict 
these values in the design and testing process. 

The simulation of the nozzle exhaust using combustion products in a conventional wind 
tunnel is quite difficult, and is often performed with a substitute gas. Wind tunnel tests are being 
conducted at NASA Langley Research Center using several different simulant gases, including 
air and mixtures of Argon and Freon. It has been shown that the proper mixture of Argon and 
Freon can simulate the specific-heat-ratio (7) of the hot scramjet exhaust products, in addition 
to matching the correct values of Mach numbers, pressure ratios, and temperature ratios [1—3] . 
This study is concerned with the computational analysis of the scramjet afterbody flowfield using 
air as the simulant gas (see Appendix A for a description of the simulant gas concept). This study 
was conducted in parallel with a computational analysis of the same scramjet nozzle- afterbody 
model using an Argon-Freon mixture as the simulant gas [15-17]. 

A wind tunnel model of a single module scramjet nozzle-afterbody configuration was 
constructed for testing in the NASA Langley 20-inch Mach 6 tunnel. The maximum Reynolds 
number of this tunnel is approximately 23xl0 6 /meter with stagnation pressure and temperature 
of 2.5 MPa and 480K, respectively. In Fig. 2.1, a two-dimensional cross sectional drawing of 
the model is shown. The simulant gas mixture is fed into a high pressure plenum chamber via 
a mounting strut. The gas in this plenum chamber is expanded through a converging-diverging 
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supersonic nozzle to approximately Mach 1.7 at the combustor exit plane (location x 3 in the 
figure), where it is further expanded over the nozzle-afterbody section of the model. This 
supersonic exhaust flow also encounters a hypersonic (Mach 6) freestream air flow, through 
which mixing occurs in a free shear layer containing additional expansions and shock waves. 
Instrumented pressure taps are located along the external nozzle-afterbody surface to measure 
the static pressures over the ramp. In Fig. 2.2 an isometric drawing of the model is shown. A 
removable tapered flow fence is shown which, when in place, is used to simulate a quasi two- 
dimensional flow. When this fence is removed, the nozzle flow also mixes with the hypersonic 
freestream in the lateral direction through a spanwise expansion, causing the flow to become 
fully three-dimensional. 

The design and analysis processes for this type of a nozzle-afterbody section is very complex 
due to the fact that many additional parameters must be considered, in addition to those which 
must be accounted for in conventional nozzles. This particular nozzle is highly asymmetric, 
and consists of an internal and an external portion. The forces and moments generated by most 
conventional nozzles can be determined by analyzing the flow up to the nozzle exit plane only. 
In this particular case, the analysis must extend further downstream due to the fact that the 
lower aft portion of the aircraft forms the external portion of the nozzle. The flow over this 
afterbody region will have a dramatic effect on the thrust vector and pitching moment generated 
by the engine module. 

The location and structure of the shear layer also becomes a concern. The internal and 
external flow, which are initially separated by the nozzle cowl, join together and mix through 
a shear layer. In addition to the usual fluid and thermodynamic properties which vary across 
the shear layer, there is a large pressure gradient between the two flows, which causes the shear 
layer to bend inward or outward depending on the sign of the pressure gradient between the 
internal and freestream flows. 

An additional flow phenomenon encountered in this study is flow separation. Flow separation 
occurs when the momentum of the flow within a boundary layer region becomes too small to 
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overcome a positive (adverse) pressure gradient (i.e. > 0). This results in the flow breaking 

away from the boundary, causing a backflow or reversed flow region. The actual point of two- 
dimensional separation can be loosely defined as the point where the velocity gradient near the 
surface in the direction normal to the boundary (§^)w becomes zero. 

In the general case of supersonic nozzle flow, there are several exhaust conditions which 
may exist. Among these are overexpanded (off-design) flow at the nozzle exit plane, and 
underexpanded flow at the nozzle exit plane. These two exhaust conditions are discussed below. 
In the internal region of the nozzle, the upstream flow at the nozzle throat feeds downstream 
towards the afterbody surface and undergoes a centered expansion as it reaches the lower comer 
of the ramp. Another expansion occurs as the flow encounters the upper corner of the nozzle 
cowl. The wall boundary layers along the cowl surfaces become an expanding shear layer as 
the flow moves downstream past the nozzle exit plane. 

In the case of overexpanded flow for the asymmetric nozzle of this study, as the flow reaches 
the nozzle exit plane, the static pressure is lower than that of the external freestream pressure, 
which results in a shock emanating from the tip of the nozzle cowl. A two-dimensional sketch 
of this type of flow is shown in Fig. 2.3. This type of flow would be highly undesirable in an 
actual flight condition since the shock wave impinges on the afterbody surface of the airframe. 
This could have an adverse effect on the stability and trim of the aircraft. 

An underexpanded flow through the same asymmetric nozzle, results in the internal nozzle 
flow continuing to expand rapidly down the afterbody surface. In this case, the upstream internal 
nozzle flow goes through the same two centered expansions at the lower ramp comer and the- 
upper cowl comer, but when the flow reaches the nozzle exit plane, instead of encountering 
a shock, it goes through the additional expansion to match the freestream conditions. At the 
tip of the cowl the interaction between the expanding jet and the external flow produces a lip 
shock and a contact discontinuity. The two-dimensional structure of this flow is shown in the 
sketch of Fig. 2.4. 
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If the internal nozzle flow is allowed to mix with the freestream flow laterally, the flow 
becomes three-dimensional, and the same type of shocks and expansions can occur in the 
spanwise direction. In this study, the two-dimensional overexpanded and underexpanded nozzle 
flows are examined computationally. The two-dimensional asymmetric nozzle is then extended 
in the lateral direction, from the sidewall reflection plate up to the spanwise plane of symmetry 
shown in Fig. 2.2. For this case, the flow structure is nearly the same as that of the two 
dimensional case (since the nozzle-afterbody region is spanwise symmetric), with the addition 
of a viscous boundary layer along the sidewall reflection plate. The flow in this case will remain 
essentially two-dimensional (i.e. there is no lateral expansion or compressions which would send 
expansion waves or shocks across the afterbody surface). 
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Chapter 3 

FORMULATION AND GOVERNING EQUATIONS 

The three-dimensional governing equations are given in detail in Sec. 3.1. The details of 
the Baldwin-Lomax turbulence model are discussed in Sec. 3.2, and the modifications to the 
model for this particular problem are given in Sec. 3.3.. 


3.1 Governing Equations 

The equations used to describe the three-dimensional viscous flow of a compressible fluid 
are the time averaged Navier-Stokes equations. These equations in terms of mass averaged 
conserved variables can be written in vector cartesian coordinate form as: 
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where Stokes’ hypothesis (A = -2/3 p) has been assumed. Heat flux terms are formed from 
Fourier’s heat conduction law. 
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The molecular viscousity, p\, is assumed to be a function of temperature only, and is computed 


using Sutherland’s formula, 



5 T 0 + S 
T + S 


(3.5) 


where S is the Sutherland constant and po and To are the reference viscosity and temperature, 
respectively. The effects of turbulence are modeled by adding an eddy viscosity, p t . This eddy 
viscosity is computed from an algebraic turbulence model and is discussed in detail in Sec. 3.2. 
The total energy and internal energy are given by 




+ ^p(u 2 + V 2 + w 2 ) 


and e = C V T 


(3.6) 


The equation of state used to close this system of equations is given by the simple ideal 
gas law. 
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(3.7) 
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For convenience, the governing equations are recast in non-dimensional form. The dependent 
variables have been non-dimensionalized as follows: 
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where the non-dimensionalized variables are denoted by an asterisk. The oo subscript indicates 
a freestream value, and / is a reference length. For convenience sake, the asterisks on the 
non-dimensionalized variables will be dropped henceforth. 

These governing equations have been written in physical cartesian coordinates (x,y,z). They 
can however, be transformed to generalized curvilinear (body fitted) coordinate space (£,t/, 0- 
The curvilinear coordinates are related to the cartesian coordinates by: 
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The transformed form of Eq. (2.1) can be written as 
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For further details on this transformation see Ref. 32. The inviscid and viscous flux vectors in 
the generalized curvilinear coordinate system are defined as 
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where U, V, and W are contravariant velocities defined as 
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These equations can be greatly simplified for a two-dimensional flow by dropping the H and 
H v vectors, and neglecting all derivatives in the ( direction. 
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3.2 Turbulence Model 


In order to properly model high Reynolds number flows, some method must be applied to 
predict the effects of turbulence within the fluid. The time averaged Navier-Stokes equations, 
given in Sec. 3.1, will only resolve the mean characteristics of the flow; they do not contain 
enough information to completely resolve the turbulent structure of the flow field. Therefore, 
a turbulence model is added to the set of governing equations to account for the effects of 
turbulence. 

To this date, many different turbulence models have been formulated, none of which has 
been deemed universal to all types of flow. The biggest problem in modeling turbulence lies 
in the fact that many different length scales exist within a turbulent flow field. The two layer 
algebraic model of Baldwin and Lomax [21] is a widely accepted model due to its reasonable 
approximations of turbulent effects, and its ease of implementation in to a finite difference type of 
algorithm. Therefore, this is the model that has been chosen for this study. Several modifications 
have been made to this model to account for special conditions which exist in the flows of this 
study; these modifications are discussed in Sec. 3.3. 

The Baldwin- Lomax model simulates the effects of turbulence by using an artificial eddy 
(turbulent) viscosity. The Reynolds stresses which arise from the time averaged Navier-Stokes 
equations are assumed to be proportional to the laminar stress tensor, where the coefficient of 
proportionality is defined as the eddy viscosity (/i t ). The Baldwin-Lomax model is a two layer 
model, and defines the eddy viscosity as: 

(/^)» nner , y<y 

fit = { (3.16) 

(/dottier , y>ycTosaovtT 

where y is the normal distance from the wall, and y CT ossover is the smallest value of y at which 
the values from the inner and outer formulas are equal. 

In the inner region, the Prandtl-VanDriest mixing length formulation given as 

Minner = P P M (3J7) 
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where L>l is the magnitude of the local vorticity vector, and 


here. 


/=fcy[l-ezp(jf + / J 4 + )] 


(3.18) 


-f y PwTw 

y = — 


{subscript w denotes wall values) 


(3.19) 


In the outer region 


(Pouter == K C cp pF wake F k i eb (y) 


with 


Fuynhe 


ke = min < 


ymax F max 

C wk ymax w diff / Fmax 
where F max is the maximum value of the function F(y) defined as: 




(3.20) 


(3.21) 


F(y) = y\u\ [l - exp(-y + /A + )\ (3.22) 

and y ma x is the value of y at which F max occurs. The function F k i e b(y) is the Klebanoff 
intermittency factor, and is defined as: 

Fkleb(y) ■ 

Udifr is the difference between the maximum and minimum total velocity at a fixed x (or £) 
station. 


1 + 5 . 5 ^^^) 

\ Vmax J 


-1 


(3.23) 


u diff = (\/u 2 + U 2 + U? 2 ) ~ (y/v? + U 2 + W 2 } 

V / max V / mir, 


(3.24) 


The second term in Udiff is set to zero (except in wakes). 

The constants which appear in Eqs. (2.21-27) are given the following values (from [21]): 

4+ = 26 


C cp = 1.6 
Ckleb = 0.3 
k = 0.4 


I< = 0.0168 
C w k = 0.25 
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One major advantage of this model lies in the fact that the location of the boundary layer edge 
does not need to be computed, since the length scales are based on the distribution of vorticity. 

3.3 Turbulence Model Modifications 

The major difficulty encountered in applying the Baldwin-Lomax turbulence model to flows 
where regions of separation exist lies in the evaluation of a suitable length scale, y m ax> and in 
turn determining (/< t ) 0 uter in the separated region. In Fig. 3.1a, a general F(y) curve is shown 
(Eq. (3.22)). If a region of separation or an overlying vortex exists, the F(y) curve may switch 
to that of Fig. 3.1b. In addition to the local peak in F(y) at y=a from the attached flow region, 
the separated flow region causes another peak to occur at y=b. The choice of this peak at y=b 
due to the separated flow results in an erroneous value of F wa k e , and in turn, a value of (^t)outer 
which is much too high. Thus, in general, the computed value of the eddy viscosity in the 
separated region will be too large, causing the details of the flow in this region to be washed out 
or distorted. In order to alleviate this problem, the Baldwin-Lomax turbulence model has been 
modified using the strategy of Degani and Schiff [22]. 

Degani and Schiff proposed cutting the F(y) curve off after the first peak is reached, so as 
not to pick the higher value of the second peak. In this method, the code searches for the first 
peak in F(y) outward from the wall to the free stream along each computational coordinate. A 
peak is considered to be found when the value of F(y) drops to 90% of the local maximum value. 
The value of 90% is somewhat arbitrary, and is chosen based on the type of flow being solved. 
In this case, the peaks in F(y) are spaced far enough apart such that the logic described will 
pick up the first peak. However, if additional separation occurs in the vicinity of the primary 
separation (e.g. crossflow separation in a three-dimensional flow), the second peak might be 
chosen, since F(y) does not drop to 90% of the local maximum after the first peak. This may 
be avoided by specifying a cutoff distance in terms of y ma x from the previous local y max , i.e. 
ymax = Cy ma x(previous)> where C is a constant chosen equal to 1.5. 
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The second modification made to the standard Baldwin-Lomax model is the inclusion of 
multiple wall effects, for computing the value of fi t near comers [23]. In this case, the eddy 
viscosity was computed as a singular function from each wall, and then an effective eddy 
viscosity was computed by taking an inverse average of the values computed from each wall. 
For example, on the comer bounded by the lower ramp surface and the viscous side wall, the 
eddy viscosity was computed as: 


fk 


(#•)<, + (fe)„ 

v/(!/ + )r, 2 + 


(3.25) 


Here the lr and sw subscripts denote lower ramp and side wall values respectively. The greater 
influence on the eddy viscosity is thus determined by the wall with the lower y + value. 

The final modification made to the model is done so to account for the influence of the 
turbulence generated by the cowl which propagates into the shear layer. This is approximated 
by using a relaxation eddy viscosity model to represent the different length scales in the problem. 
Following the work of Waskiewicz, et. al., [24], the eddy viscosity in the wake is computed as: 


(k = (fit) c + + (MeH 1 “ exp(-x c //3S c )\ (3.26) 

where (/x t ) c is the calculated eddy viscosity at the tip of the cowl, and (// t ) v is the calculated 
“outer” eddy viscosity value based on the local vorticity value in the wake. The distance between 
these two stations is denoted by x c , and S c is the instantaneous boundary layer thickness at the 
upper cowl tip. The parameter is a relaxation length scale given a value of 10. 
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Chapter 4 

SOLUTION ALGORITHM 


The solution algorithms for both the two-dimensional and the three-dimensional governing 
equations are based on implicit finite- volume methods. In Sec. 4.1 the basic finite volume 
methodology is described. The time integration of the set of discretized equations is given in 
Sec. 4.2, and the implementation of the physical boundary conditions is discussed in Sec. 4.3. 
The computational grid generation methodology is discussed in Sec. 4.4. 

4.1 Finite Volume Discretization 

Both the two-dimensional and the three-dimensional sets of governing equations are solved 
using implicit finite volume schemes. In this method the integral formulation of the governing 
equations in conservation form is discretized directly in the physical space. This direct method 
of discretization ensures the conservation of mass, momentum, and energy at the discretized 
levels, even in the presence of strong discontinuities such as shocks. Equation (3.10) can be 
expressed in integral formulation as 

ITtJjJ ® dv + J J F hds = 0 (4J) 

_ V 3 

where the vector F is defined as 

F=(e- E v ) i+(F- F v )j +(g~ G v ) k (4.2) 

and nis a unit normal vector pointing outward from the surface S bounding the volume V, 

h = n x i + n y j + n z k. (4.3) 

The restriction of these equations to the two-dimensional equations is straightforward. 
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The flowfield is discretized into a set of ordered cells (volumes), and Eq. (4.1) is applied 
to each volume. The semi-discrete form of Eq. (4.1) at a cell whose location is i, j, k can 
be written as 


dQ 

dt 
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= 0 


(4.4) 


In the computational domain AC, A? 7 , and AC are arbitrary, and for convenience taken as 


AC — Ci+l/2 — Ci— 1/2 — 1 


(4.5a) 


Ar/ = Vj+l/2 ~ Vj- 1/2 = 1 


(4.5b) 


AC = Cjt+1/2 - Cifc-1/2 — 1 


(4.5c) 


4.2 Time Integration of the Discretized Equations 

The system of governing equations described in Sec. 4.1 is solved using an implicit, finite- 
volume, upwind, spatially factored scheme [25, 26]. This scheme employs upwind differencing 
for the convective and pressure terms, and central differencing for the diffusive terms. The 
upwind differences are constructed using a Roe flux-difference splitting method [27, 28]. 

Since the flows considered in this study are primarily steady-state in nature, a psuedo- 
time marching approach is used where the governing equations are treated as time dependent 
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Navier-Stokes equations. These equations are hyperbolic in time and allow for time marching or 
propagation type solutions.Equation (4.1) must first be linearized in time. The non-linear terms 
occur due to the fact that the flux vectors E, F, and G at the n+1 time level are functions of 
the vector of conserved variables Q. This linearization is accomplished by expanding the flux 
vectors at (n+1) in a Taylor series about (n) as. 


(4.6a) 


= kn + \ ad ) (^” +1 ~ A< + A<)2 


A E n = A n AQAt + ti(Aty 
where A is the flux Jacobian matrix |£. Similarly 


(4.6b) 


A F n = B n AQ n At + 0(A if 


A G n = C n AQ n At + tf(A t) 2 (4.8) 

where B and C are the flux Jacobian matricies ^and respectively. The solution is marched 

oQ dQ 

in time by applying an implicit Euler integration scheme [26] to Eq. (4.4). Substituting Eqs. 
(4.6b), (4.7), and (4.8) into Eq. (4.4), and simplifying. 


T i O rv 

JS + ZpSo = -*«?") 


where A Q = Q n+1 — Q n , R n is the steady-state residual evaluated at time level (n), and 

[a^-i.) a(F-/„) a(6 -<?„)' 

R= Ti + 5, + a< ( 

Substitution of Eqs. (4.6-8) into Eq. (4.9) yields: 


(4.10) 


J At + 6 t dQ 


d(P-F v ) d(G-G v ) d( H-H v \ 
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When Eq. (4.11) is written at discrete points in the computational domain, it results in a large 
banded matrix, which is expensive to solve computationally. One method of circumventing 
this problem is to use an approximately factored form of the equation. A number of different 
methods of factorization are possible for Eq. (4.11), one of which is the three-factor block- 
tridiagonal scheme [29] in which the implicit operator is split spatially such that each factored 
set contains flux Jacobians for its respective direction only. Applying this factorization to Eq. 

(4.11), one obtains 


I 

J At 

I 

J At 



I d(6-G v ) 

L 1 C_ 

J At v dQ 

A Q n = Res 11 


(4.12) 


This method results in a set of block-tridiagonal matricies which can be solved in three 
computational sweeps. In the first sweep, the ^-direction differences are treated as implicit, 


I 

JAt 



(4.13) 


where 



(4.14) 


Equation (4.13) is solved for A Q* . Equation (4.14) is then rearranged for the ^-direction 
implicit sweep, 


where 



(4.15a) 


(4.15b) 
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Equation (4.15) is solved for AQ**. Then Eq. (4.15b) is rearranged as follows, in order to 


solve for AQ n , 


7Ei + s < 


d(d-G v ) 

OQ 


A Q n =[jAt AQ ' 


The vector of conserved variables is then updated to the next time level as: 

Q n+ 1 = Q n + A Q n 


(4.16) 


(4.17) 


The upwind differences for Eqs. (4.13—16) are constructed using the flux-difference splitting 
method of Roe [26, 27]. In this method, the spatial derivatives are written conservatively as a 
balance across a volume element. For example in the £ — direction the F flux difference at point 
i (holding the j and k indices constant) can be found as: 

( di' \ _ £<±1/2 - ft-l/2 (4.18) 

where F, +1 / 2 and F i _ l / 2 are cell interface values, and can be constructed as: 

^.+i/2 = \ [^(< 3 ,"+ 1/2) + P ( Qt + 1/2) - l^li-1-1/2 (^i+1/2 _ ^i+1/2)] i+1/2 (4,19a) 


l^li+l/2 ~ '^>+1/2 li-t-l/2^i+l/2 


(4.19b) 


Here, A is the flux Jacobian matrix of F, A is the matrix of eigenvalues of A, and T and T^ 1 
are the right and left eigenvectors respectively. The overbars on the above variables denote 
Roe-averaged variables. Q + and Q - denote conserved state variables on cell interfaces, which 
are obtained from the primitive variables as: 


Q7+1/2 = Q* + j^K 1 - v* + (1 + K z) A t\Qi 


(4.20a) 


Qt + 1/2 = Qi + 1 - ^[(1 + h ) v * +(1 - 1 


(4.20b) 
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where 


A ( Qi = Qi+i - Qi and Vc Qi = Qi ~ Qi- 1 (4.21) 

The order of differencing is determined by <^; for first order differencing <j)^ is set zero, and for 
higher order differencing is set to one. The type of differencing is controlled by the parameter 
«. « = — 1 corresponds to second order fully upwind differencing, k = 0 corresponds to second 
order upwind biased, k = +1 to second order central differencing, and k - 1/3 to third order 
upwind biased. 

In order to accelerate the convergence of the three-dimensional computations, two techniques, 
known as mesh sequencing and multigridding are employed [26]. A sequence of grids are formed 
over the computational domain by deleting every other grid line on the next finer grid. In mesh 
sequencing, the computations are started by iterating on the coarsest level grid, until the flow 
solution develops. Then, that solution is interpolated onto a finer level grid, where more iterations 
are performed, with an increased level of spatial resolution. This process is repeated until the 
finest level grid is reached ( in this study, three levels are used). The rationale behind this mesh 
sequencing strategy lies in the fact that in the initial computations, when spatial resolution is not 
of vital importance, the calculations are performed on very coarse grids which helps to accelerate 
the development of the flowfield due to the large spatial steps. 

Once the flow solution is developed on the finest level grid, the multigrid strategy begins. 
The multigridding technique uses the same coarser level grids which are created in the mesh 
sequencing progression, however instead of iterating on a coarser level and then moving up to 
the finer levels, the iterations are performed at the finest level first, and then cycled through to 
the coarser levels and back up again. There are several different cycling strategies which may 
be used, several of which are discussed in detail in [26], The idea behind multigridding lies in 
the fact that the low frequency error components which occur on a given fine grid level can be 
reduced to high frequency errors by decreasing the grid resolution (i.e. using a coarser level 
grid). These high frequency errors on the coarser level grid are then reduced by applying the 
solution algorithm, which is much more capable of resolving high frequency errors. This process 


25 


increases the overall robustness of the computational scheme, which allows for larger time steps 
to be used, and hence faster convergence rates. 

4.3 Boundary Conditions 

The initial and boundary conditions play an important role in the accuracy of the solution for 
a physical flow. Therefore, special attention must be paid to the formulation of these conditions. 

The upstream boundary conditions for the nozzle and external regions are generated using 
the two-dimensional boundary layer profile program of [30]. In the internal nozzle region of the 
experimental model, the flow is expanded from a high pressure plenum through a convergent- 
divergent nozzle to a specified Mach number. However, there is no experimental data for 
the boundary layer thickness at this location (which corresponds to the upstream plane of 
the computational domain), so the computational boundary layer profiles for this location are 
generated in a somewhat arbitrary manner. The sum of the profile thicknesses from the lower 
ramp surface and the upper cowl surface inside the nozzle is assumed to be approximately 
13 percent of the nozzle height. These boundary layer profiles are generated using the two- 
dimensional boundary layer program (assuming fully turbulent flow), and then mapped onto 
the computational grid’s initial plane. The external flow profile is generated using the same 
two-dimensional boundary layer program by assuming a viscous turbulent flow over an 18” flat 
plate, which corresponds to the length of the upper cowl surface on the experimental model 
(Fig. 2.2). This results in a boundary layer thickness of approximately 0.5 inches. Were this 
code to be used as an actual design code for a nozzle-afterbody, it is likely that more detailed 
data would exist for the upstream profiles (either from more extensive experiments or from other 
computational results). 

The upstream profiles for the three-dimensional calculations are generated using the same 
boundary layer profiles of the two-dimensional cases. These two-dimensional profiles are 
“stacked” along grid planes in the spanwise direction to form upstream profile planes. The 
presence of the third viscous sidewall (in both internal and external regions) causes a third 


boundary layer to develop, along with a very weak leading edge shock. In the actual corner flow 
regions, vorticies may exist; however in the present calculations these vorticies are not modeled. 
A somewhat more realistic method of simulating the upstream three-dimensional boundary layer 
profiles was used in [17] , although once again, as in the case of the two-dimensional calculations, 
there is no experimental data regarding the boundary layer thicknesses at these locations, so the 
profiles are somewhat arbitrary. 

On regions bounded by solid surfaces, no-slip conditions are applied such that all velocities 
vanish at the boundary. These solid surfaces are also considered to be adiabatic. The pressure 
on a solid boundary is found from a zeroth -order-extrapolation from the interior point value. 
The density is then computed from the state equation. 

u = v = w - 0 , = 0 , ir- = 0 (4.21) 


The conditions at the outflow boundaries are set depending on the nature of the flow, be 
it subsonic or supersonic. If the flow is supersonic, all of the characteristic flow lines point 
from the inside of the computational domain outward, so that all of the variables at the outflow 
boundary may be determined from an extrapolation of the interior flow variables. If the flow 
is subsonic (e.g. within boundary layers) one of the characteristic lines is incoming (i.e. from 
outside of the computational domain inwards). Therefore, one of the boundary conditions must 
be specified analytically. In this case the freestream pressure is used. The remaining variables 
are extrapolated from the interior domain as before. 

The boundary conditions on the plane of symmetry (three-dimensional case) are defined by 
extrapolation of all of the variables from the interior domain with the exception of the velocity 
component which lies perpendicular to the plane of symmetry. This normal component is set 
equal to the negative value of its interior point, thus forming a symmetry plane. 

In this study, the two-dimensional computations are performed on a single body fitted grid. 
In order to allow for the application of the standard block tri-diagonal inversions from one end of 
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the computational zone to the other in each direction (r/= 0 to r/ max for //-implicit and £-sweep; 
£=0 to £ max for implicit and tj sweep), some cells must coincide with the region bounded by 
the cowl. These cells are computed as if fluid were in this region, but they are then blanked out. 
This “blanking” is accomplished by (a) setting the off-diagonal blocks (which are 4x4 matrices) 
in the coefficient matrix to zero, (b) setting the diagonal elements of the diagonal blocks to 
unity, (c) setting the residual vector components to zero, (d) the Q values of the cells which 
neighbor fringe cells (cells which lie adjacent to “blanked cells” ) are set equal to the Q values 
of the fringe cells, in order to avoid incorrect flux computation at the cowl boundary. Thus, the 
specified wall boundary conditions for these cells are automatically preserved. 

The three-dimensional computations are performed using a block structured approach. In 
this method, the physical space is divided into separate, distinct zones (in this case 4 zones have 
been used), where each block is treated as a separate computational domain. On blocks which 
lie next to each other, values of conserved variables are transferred across the interfaces along 
contingent cells. For example, the downstream boundary of one block becomes the upstream 
boundary of another. Conserved variables from two cells neighboring the boundary of one block 
are passed to the adjacent block as ghost cells, thus retaining the second order spatial accuracy 
of the scheme. The main drawback to this method lies in the fact that the lengths of the data 
vectors are restricted by the lengths of the blocks, which results in slower computational rates. 

4.4 Grid Generation and Adaptive Grids 

The physical domain over which the flow field is being sought needs to be discretized 
into a finite set of points, called the grid, to which the solution algorithms can be applied. 
Several methods exist for generating the computational grid, among them are algebraic methods, 
partial differential equation solvers (elliptic or hyperbolic), and complex variable transformation 
techniques. The proper choice of a grid is vital in generating accurate computational results 
for any flowfield. A grid which is not well suited to the flow can lead to instabilities in the 
computations or even a lack of convergence of the solution [33]. Grid density should be very high 
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in regions where large gradients will exist (e.g. boundary layers, shocks, shear layers); however, 
the number of grid points used to discretize the flow field is limited by computer memory and 
speed, so a judicious choice of grid spacing is crucial. In this study, the domain is initially 
discretized using an algebraic grid generation scheme due to the ease of implementation and the 
fact that the physical boundaries of the flow domain are easily described by explicit algebraic 
functions. All of the clustering is based on an exponential stretching/clustering function given 
in [33]. In Fig. 4.1 the initial two-dimensional grid is shown. The grid point clustering on this 
grid is generated based on no a priori knowledge of the computational results. However, the 
basic flow physics described in Chap. 2 are assumed, and the clustering of grid points is seen 
to reflect this hypothesized flow. Along the ramp surface and around the cowl walls the grid is 
very dense in the normal direction in order to properly resolve the boundary layers. The grid is 
also clustered longitudinally near the comers on the ramp surface and on the inner cowl surface 
since it is expected that centered expansions will occur at these locations. 

In addition to the computations on the fixed grid, flow adaptive grid calculations are 
performed for the two-dimensional cases. Once the solution develops on the fixed grid, the 
clustered and stretched regions of the grid are redefined based on the locations of large gradients 
in the flow. The method used to adapt the grid to these flow gradients is described in detail in 
Appendix B . The flowfield is then recomputed on the adapted grid. 

The three-dimensional grid is generated using the same clustering strategy as that of the 
two-dimensional fixed grid, with additional clustering in the spanwise direction near the viscous 
reflection plate. As mentioned before, the three-dimensional solver utilizes the block structured 
approach, so actually, four grids are generated (one for each zone). The three-dimensional grid 
is shown from several perspectives in Fig. 4.2. A three-dimensional cutaway view of the entire 
grid is shown in Fig. 4.2a. The spanwise clustering can be seen on the upstream vertical plane. 
The inner nozzle region up to the cowl tip is shown in the enlarged view of Fig. 4.2b, and the 
spanwise symmetry plane, which is essentially the same as the two-dimensional grid, is shown 
in Fig. 4.2c. 
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The number of cells in each direction and the total number of cells of each grid are tabulated 
in Tabel 4.1. The description of each case (Case 1, 2, 3) are explained in the following chapter. 


Chapter 5 

RESULTS AND DISCUSSION 

As mentioned previously, several different flow cases have been solved computationally. 
These cases will be referred to as: 

Case 1: Two-Dimensional Overexpanded Flow (off-design) 

Case 2: Two-Dimensional Nozzle Flow (design) 

Case 3: Three-Dimensional Nozzle Flow (design) 

The first case was initially intended as a trial run, to be used for code set up and initial 
validation of the grid. Although no experimental data exists for this case to compare the 
numerical solution, the results are being shown here due to the interesting nature of the flow. 
The second case was intended to be used as a benchmark of the two-dimensional code in that 
the computational results could be compared directly with the experimentally obtained results. 
The third case, which is simply an extension of the second case from two-dimensions to three 
dimensions, is also to be used in the validation process of the three-dimensional code set up. A 
summary of the upstream flow parameters for each case is given in Table 5.1. 

Since the flows being solved are assumed to be steady state in nature, a local time stepping 
strategy is used. Each computational cell is advanced in time based on its maximum allowable At 
value. In order to stabilize the initial numerical transients, the time steps are gradually increased 
by increasing the Courant number from 0.01 to about 4.0. Table 5.2 shows the number of 
iterations and total CPU time required to obtain converged solutions. In the study, convergence 
is defined as the point where the logarithm of the residual (right-hand- side of Eq. (4.11)) is 
decreased by approximately four orders of magnitude. 

The two-dimensional cases were run on the VPS-32 (an enhanced version of the Cyber 205) 
computer at NASA Langley Research Center. The three-dimensional runs were performed on 
the CRAY-2 computer, also located at NASA Langley Research Center, and used the mesh- 
sequencing/multigrid algorithm [26] described briefly in Chap. 4 to accelerate the convergence 
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rate. It should be mentioned here that the computational rates achieved with these codes are 
very good in comparison with those required for a similar analysis using an explicit, finite 
difference, MacCormack scheme [15-17]. The multiple species analyses, which are reported in 
[15-17], are performed by assuming a mixture of four gases (Argon, Freon-13, Nitrogen, and 
Oxygen). This assumption leads to additional species continuity equations which are coupled 
to the global continuity, momentum, and energy equations. The two-dimensional, explicit, 
multispecies analysis code uses more than twice as much CPU time as does the implicit single 
gas analysis code used in this study. The additional CPU time required by the three-dimensional 
multispecies code is even more drastic. Almost ten times more CPU time on the CRAY-2 
computer is needed to obtain converged solutions. The nozzle expansion rates for air are slightly 
higher than those for the multispecies stimulant gases, and thus the pressure distributions over 
the afterbody surface differs somewhat. However, for preliminary configuration analysis, the 
single gas, with constant specific-heat-ratio assumption seems to be a much more economical 
alternative to the multiple species studies. 

In Figs. 5.1-2, the residual vs. iteration history is shown for Cases 1-2, respectively. In 
both r ’-'ses 1 and 2 thr —sidual is reduced by approximately three orders of magnitude after 
about 7000 iterations on the fixed grid. Further iterations do not result in any additional decrease 
of the residual at this point. After 10,000 iterations, the grid is adapted to the flowfield, using the 
grid adaption scheme described in Appendix B. The error which results from the interpolation 
of the fixed grid solution onto the adapted grid causes the sharp increase in the residual values 
at this point. After an additional 6000 iterations, the residual is reduced by approximately four 
orders of magnitude, at which point it begins to level out again. Additional iterations at this 
point do not decrease the residual any further. 

The residual history of Case 3 is shown in Fig. 5.3. The first 5000 iterations are performed 
on the coarsest level grid, at which point the residual is reduced by about one order of magnitude. 
The spikes which occur up to this point are the result of restarting the solution from a previous 
iteration level. This restarting strategy is required due to a CPU time constraint. After 5000 
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iterations on the coarsest level, the solution is “prolonged,” or interpolated, onto the next finer 
level grid. An additional 1500 iterations are performed on this level, at which point the solution 
is “prolonged” to the finest level grid. The sharp increases in the residual values after 5000 and 
6500 iterations occur due to the errors introduced by the interpolation from the coarser level 
grid to the finer level grids. The multigrid cycling strategy is invoked at this point, and the 
residual is seen to drop dramatically to the point where it is reduced by approximately four 
orders of magnitude. 


5.1 Case 1 

The results of this case 1, which depict the overexpanded nozzle flow, are presented in Figs. 
5.4-8. This case was initially set up and run on the fixed grid ( Fig. 4.1). The specific-heat-ratio 
( 7 ) was assumed to be constant and equal to 1.214, which corresponds to the specific-heat-ratio 
of the 72% Freon 28% Argon (by mass) mixture used as a simulant gas. The Mach number and 
static pressure contours obtained from solutions on the fixed grid are shown in Figs. 5.4 and 5.5, 
respectively. The internal nozzle flow goes through a centered expansion at the lower comer of 
the throat where the 20° afterbody ramp begins. This expansion fan strikes the lower cowl surface 
just upstream of the lower cowl corner where the flow undergoes a second centered expansion. 
At the point where the nozzle flow reaches the exit plane it is overexpanded with respect to 
the freestream external flow. This overexpansion results in a shock, with a pressure ratio of 
approximately 2.7, emanating from the tip of the cowl. This shock impinges on the afterbody 
surface about 1/3 of the distance down the length of the ramp, resulting in a weaker reflected 
shock. The shock-boundary layer interaction on the ramp also results in a separation bubble, or 
a region of reversed flow. The interaction between the internal and external flow occurs through 
an expansion and a shear layer, which also emanate from the cowl tip. The top expansion 
fan points upwards with an included angle of approximately 55°. The shear layer is pointed 
downward, following the same line as the afterbody surface, and eventually interacts with the 
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reflected shock as it reaches the downstream portion of the ramp. The pressure ratio of the fully 
expanded downstream location to that of the upstream nozzle conditions is approximately 0.3. 

The separation region on the ramp is further evidenced by the plot of velocity vectors with 
the streamlines superimposed (Fig. 5.6). In addition to the separation on the ramp caused by 
the shock-boundary layer interaction, a small region of reversed flow is observed right on the 
inside of the cowl tip. This is due to the rapid expansion of the external flow around the cowl 
tip, and the lower momentum of the boundary layer flow within the nozzle along the cowl wall. 

Once the solution is obtained on the fixed grid, the adaptive gridding algorithm, described in 
Appendix B, is invoked. The grid is “adapted” to the specific flow field. Various parameters and 
combinations thereof may be used as weighting functions in the grid adaptation. In this case, a 
combination of the density ( p ) and the total velocity (V=u 2 +v 2 ) is used. In regions where large 
gradients of p and/or V occur, the grid spacing will become very dense. This is shown in Fig. 
5.7 where the adapted grid is seen to closely mimic the actual flow structure. The grid becomes 
very clustered where shocks and expansions occur, and somewhat sparse where there are no 
large gradients in velocity and density. One of the major problems encountered in applying 
die adaptive grid scheme to this particular flow geometry is the fact that the cowl boundaries 
protrude into the computational domain. If the adaption scheme was applied directly to the grid 
as a whole, the geometry of the cowl would be changed drastically, since very large gradients of 
density occur between the internal and external walls of the cowl. To alleviate this problem, the 
adaption process is split into two regions. In the first, the region from the upper cowl surface 
to the outflow boundary (y-direction) is adapted. The second region, which consists of all of 
the remaining area (i.e. everything that lies below the external boundary of the cowl surface) 
is then adapted. This method does place some restriction on the adaption of the lines dividing 
these two regions, however the flow gradients from the fixed grid solution do not appear to be 
excessively large here, so this is a minor inconvenience. 
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The Mach contours obtained from the solution on the adapted grid are shown in Fig. 5.8. 
The resolution of the solution near the shock regions is greatly enhanced over that obtained from 
the solution on the fixed grid, as can be seen by comparing Figs. 5.4 and 5.8. 

5.2 Case 2 

The results of this case 2, which depict the underexpanded nozzle flow, are shown in Figs. 
5.9-15. This case was initially run using the same fixed grid as that of case 1, then rerun 
on the flow adapted grid generated using the fixed grid results. The same splitting strategy 
is used in this adaption, as was used in case 1 , to avoid changing the geometry of the cowl. 
The final adapted grid is shown in Fig. 5.9. The calculations for this case were performed 
for air ( 7 = 1 . 4 ) flowing supersonically in the nozzle region and hypersonically in the external 
region. These conditions match those of the experimental test case [ 6 ], so that the results can be 
compared directly with one another. In Fig. 5.10, the computational and experimental values of 
the ramp surface pressure coefficient (Cp) are compared. The experimental data is taken from 
the spanwise plane of symmetry, where the flow is essentially two-dimensional. The computed 
values are seen to agree quite well with the experimentally obtained values. The computational 
values of C p are plotted on both the ramp surface and the inner cowl surface, and are seen to 
decrease abruptly at the comers, where the centered expansions occur, and gradually along the 
20° ramp and 12° cowl. Additional experimental runs are currently being conducted at NASA 
Langley Research Center to generate off surface pressure data, so that the computational and 
experimental results can be compared in the shear layer region where the adaptive grid solution 
provides more resolution. 

In Figs. 5.11 and 5.12, the density and the Mach number contours are shown from the fixed 
grid results. The nozzle flow goes through the same two centered expansions at the comers 
as the flow of case 1 ; however, when the flow reaches the exit plane it is underexpanded, and 
thus continues to expand down the ramp and outward into the freestream flow. In this case the 
shear layer, which originates at the cowl tip, is deflected upwards at an angle of 7° to 15 . This 
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deflection of the inner region flow into the outer flow results in a weak external shock which 
turns the external flow upwards to match the flow direction of the shear layer. The improvement 
in the shear layer resolution obtained from using an adapted grid can be seen by comparing the 
Mach number contours from the fixed grid solution ( Fig. 5.11) and those obtained from the 
flow adapted grid solution ( Fig. 5.13). The velocity deficit extends much further downstream 
in the solution obtained on the adapted grid. This is also seen in the comparison of the velocity 
vectors plots of Figs. 5.14 and 5.15 (the fixed grid and adapted grid results, respectively). In 
this case, a slight region of reversed flow exists on top of the cowl, due to the rapidly expanding 
internal flow. This is easily seen in the velocity vector plot of the fixed grid results. 

5.3 Case 3 

The results of this case, which depict the three dimensional spanwise symmetric underex- I 

panded nozzle flow, are shown in Figs. 5.16-19.In Fig. 5.16, nondimensional ramp surface 
pressures are shown at various spanwise locations. The computational pressures tend to agree j 

with those obtained experimentally except upstream near the sharp expansion corner. This dis- S 

crepancy may be due to the lack of sufficient grid clustering near the expansion region. The 

j i 

implicit damping which is inherent to upwind schemes tends to spread the expansion region 

i 

out over several computational cells, which causes the computational pressures to lag the actual 

| 

experimental values. The non-dimensional pressure contours, which are shown in Fig. 5.17 on 

the reflection plate and the afterbody ramp, indicate that the flow is essentially two-dimensional : 

as would be expected. The flow structure is very similar to that of the flow in Case 2. The 5 

5 - s 

I § 

centered expansions occur at the lower ramp comer and the upper cowl comer in the same f 

manner as before, and the flow continues to expand as it extends down the afterbody ramp and | 

5 i 

out into the freestream. i 

1 — ” ft 

I S 

The nondimensional density contours for Case 3 are shown in Fig. 5.18. The view on the | 

left shows the density contours on the spanwise plane of symmetry. A comparison with the j 

density contours of Case 2 ( Fig. 5.11) shows almost identical results. The centered expansions \ 
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at the two comers inside the nozzle occur in the same fashion as in the flow of Case 2. The 
afterbody flow also continues to expand down the ramp and into the external freestream in 
the same manner as the two-dimensional flow, since it is not allowed (by the symmetry plane 
boundary conditions) to expand in the spanwise direction. The view on the right hand side of 
Fig. 5.18 shows the three-dimensional effects of the viscous sidewall. The density of the gas is 
higher near the viscous walls due to the increased temperature, which is a result of the adabatic 
wall boundary conditions which are imposed there. 

The Mach number contours for case three are shown in Fig. 5.19. Again, the view on 
the left shows the Mach contours on the spanwise plane of symmetry. These contours are also 
very similar to those of the Mach contours obtained in the results of Case 2. The shear layer 
is approximately the same thickness, and is seen to be deflected outwards into the freestream 
region at nearly the same angle. However, the resolution inside the shear layer region is not as 
good as the results of Case 2, and is again thought to be the result of inadequate grid refinement. 
The view on the right shows a three-dimensional perspective view of the Mach contours. Here, 
the boundary layers and shear layers show up very well. Again, the flow is seen to be essentially 
two-dimensional in all respects, except near the viscous sidewall. The flow is seen to expand 
less rapidly near the sidewall, as evidenced by the downward dip in the Mach contours in the 
external flow region close to the wall. The wall boundary layer is essentially hindering the 
expansion into the freestream. 
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Chapter 6 

CONCLUSIONS 


The viscous mixing of a supersonic nozzle-afterbody flow with an external freestream flow 
is computed successfully for several cases. A two-dimensional overexpanded (off-design) nozzle 
flow is solved on both a fixed grid and a flow adapted grid. The results of the adapted 
grid solution show much better resolution of the recompression shock and the resultant shock 
boundary layer interaction region than those obtained on the fixed grid. The two-dimensional 
underexpanded nozzle-afterbody flowfield is also solved on a fixed and a flow adapted grid. In 
this case, the shear layer resolution is slightly better in the solution obtained on the flow adapted 
grid. The computational pressure distribution on the ramp surface compares favorably with the 
experimental surface pressures obtained at the nozzle symmetry plane. 

The three-dimensional spanwise symmetric flow is solved on a multiblock, fixed grid. The 
results obtained compare well with the experimental surface pressure data and with the two- 
dimensional symmetry plane computations. This flow is essentially two-dimensional except near 
the regions bounded by the viscous reflection plate, and hence the flow structure is nearly the 
same as that of the two-dimensional case. 

The Reynolds stresses are simulated using the Baldwin-Lomax two layer algebraic model 
with modifications to account for regions of separated flow, multiple wall geometries, and 
turbulent wake flow. There are many turbulence models available which may yield more accurate 
results. However, the simplicity of the Baldwin-Lomax model, and its ease of implementation 
into finite volume type codes, such as the ones used in this study, are the reasons for this choice. 
Any errors which may be due to the Baldwin-Lomax model should most likely show up in 
the shear layer region. There is no way to verify this, since there is currently no off-surface 
experimental data available. Also, since the occurrence and location of flow separation is strongly 
dependent on the nature of the flow, be it laminar or turbulent, the accuracy of the turbulence 
model used becomes very important if flow separation is a primary concern in the analysis. 


Both the two-dimensional and the three-dimensional schemes used in this study are compu- 
tationally very efficient methods for solving the complex flows about nozzle-afterbody sections. 
The use of an implicit based scheme, as opposed to an explicit MacCormack type algorithm, 
allows for larger time step increments to be used, and hence, faster convergence rates. The 
constant specific-heat ratio (7) assumption for the internal simulant gas, leads to some deviation 
in the flowfield from that of an actual variable 7 gas of scramjet exhaust products. But the 
savings in computational time achieved by using the constant 7 assumption makes this an attrac- 
tive alternative for use in the preliminary design stages of nozzle-afterbody configurations. The 
two-dimensional CFD capabilities which have been developed and demonstrated in this study 
should be used to augment the conventional wind tunnel studies of scramjet nozzle-afterbodies. 
However, in order to include the effects of spanwise expansions, three-dimensional calculations 
must be performed for the full-span nozzle. 

Although the computational surface pressures agree well with those obtained experimentally, 
more measurements are needed to further validate these codes. A more accurate means of 
specifying the upstream boundary layer profiles is needed. This requires some experimental data 
on the boundary layer thicknesses at the nozzle throat. Also, off surface experimental flowfield 
data is needed, especially in the shear layer regions. 

Additional three-dimensional computations are currently underway to include the spanwise 
expansion of the nozzle flow into the external freestream, along with computations using the 
Freon-Argon simulant gas mixture. Also, the use of three-dimensional flow adaptive grids 
should be studied, since the results from the two-dimensional adaptive computations show better 
resolution of the flow in high gradient regions. 
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Table 4.1 Computational grid sizes 


Case 

Block 

idim 

jdim 

kdim 

Total 

1 

1 

155 

1 

132 


2 

1 

155 

1 

132 


3 

1 

33 

25 

41 


3 

2 

49 

25 

41 


3 

3 

65 

25 

41 


3 

4 

65 

25 

41 

217,300 


Table 5.1 List of upstream conditions for three cases 


Case No. 

Mach No. 

Re/m 

T t (K) 

Pt (kPa) 

7 


(int/ext) 

(int/ext) 

(int/ext 

(int/ext) 


1 

1.7 

4.92x1 0 5 

467 

172 

1.214 


1.7 

4.92x1 0 5 

467 

172 


2 

1.665 

1.26xl0 7 

475 

1166 

1.4 


6.0 

2.27x1 0 7 

478 

2520 



6.0 

2.27x1 0 7 

475 

1166 

1.4 




478 

2520 



Table 5.2 CPU time requirements for each case 


Case 

Iterations to 

Total CPU 

Speed 

Computer 


convergence 

time 

(time/iteration/gridpoint) 


1 

4500* 

6,000 sec 

34 /zsec 

VPS-32 

2 

4500* 

6,000 sec 

34 /zsec 

VPS-32 

3 

7300* 

8,400 sec 

76 /zsec 

CRAY-2 


♦including grid adaption 
**with mesh sequencing and multigriding 
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FLOWRELD ABOUT A GENERIC HYPERSONIC VEHICLE 



Fig. 1.1 Generic Hypersonic Aircraft 


NASA Langley conventional wind tunnel model 
for single-module scramjet engine 


(sectional view) 




Fig. 2.3 Sketch of the overexpanded nozzle flow 



Fig. 2.4 Sketch of the underexpanded nozzle flow 
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Figure 3.1a Behavior of F(y) in a attached boundary layer 



Figure 3.1b Behavior of F(y) in a separated boundary layer 


Fig. 3.1 Behavior of F(y) 




Figure 4.1 Two-dimensional fixed grid 
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Fig. 4.2 Three-dimensional fixed grid 
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Fig. 5.1 Case 1 Residual vs. Iteration 
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Fig. 5.2 Case 2 Residual vs. Iteration 
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Fig. 5.3 Case 3 Residual vs. Iteration 
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Figure 5.6 Case 1 streamline/velocity vectors (fixed grid) 



Figure 5.7 Case 1 adapted grid 
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Figure 5.8 Case 1 Mach contours (adapted grid) 
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Figure 5.9 Case 2 adapted grid 
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Figure 5.16 Case 3 spanwise surface pressure contours 
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Fig. 5 17 Case 3 pressure contours 
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;. 5.18b Case 3 Density contours at various crossflow 
planes 
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Appendix A 

COLD GAS SIMULATION OF HOT EXHAUST PRODUCTS 

In order to accurately simulate the nozzle exhaust flow from scramjet engine modules using 
scaled test models, several nondimensional parameters must be matched between the test flow and 
the exhaust products of the actual scramjet engine. For inviscid flow similitude, the parameters 
which must be matched include the geometry scales between the actual nozzle afterbody and 
the model, the ratio of static pressures (internal to external), the Mach numbers, and the flow 
directions at geometrically similar locations within the flow fields. If these parameters are all 
matched at the nozzle exit plane, the wave systems and Mach lines which emanate from the test 
model nozzle will match those of the actual engine. 

If a simulant gas is used, in lieu of the hot scramjet exhaust products (which are very difficult 
to use in a conventional wind tunnel), the ratio of specific heats (7) must also be matched [1] 
so that the wave systems will travel at the same speed. This is due to the fact that the speed of 
sound in a gas is directly related to 7 by the relationship a=7RT. If the pressure distribution on 
the afterbody surface of the model is to match that of the actual nozzle-afterbody, the system 
of expansion and compression waves must be similar. 

The use of simulant gases can be used to simulate the flow of hot scramjet exhaust products 
at greatly reduced temperature levels, and allow for the scaled scramjet nozzle-afterbody model 
testing in a conventional wind tunnel. 

Several different gases have been identified for use as a scramjet exhaust simulant gas [ 1 ], 
some of which are simple mixtures of Freon and Argon. In Fig. A.l the variation of 7 with 
temperature is shown for a 70 % F13B1 (Freon- 13 ) + 30 % Argon mixture (by volume), along with 
the 7 variation of a scramjet exhaust gas. As seen in this figure, the 7 values of the Freon- Argon 
mixture match those of the exhaust gas at a significantly lower temperature. 

When air is used as the simulant gas, some error is introduced into the experimental pressures 
due to the difference in specific heat ratios. The gamma value for air is approximately constant 
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at 1.4 over the temperature which is suitable for conventional wind tunnel testing. However, air 
is much less costly to use in a wind tunnel test, so it remains a good candidate for preliminary 
analysis. 

The requirements for viscous similitude are somewhat more restrictive than those of the 
inviscid flow. In addition to the aforementioned parameters, the Reynolds number, Prandtl 
number, and Schmidt number must also be matched. In this study however, the effects of surface 
heat transfer and the diffusion of the exhaust gases into the freestream are not considered, so 
it is not necessary to match the Prandtl and Schmidt numbers, and only the Reynolds number 
becomes important. 



1 1 1 1 1 I i i i 

1000 1400 1000 2200 2600 

T°K (for Scramjet Exh-gas only) 


Figure A. I Matching of specific-heat-ratios between simulant gas and exhaust prod- ucl $ 
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Appendix B 
ADAPTIVE GRIDS 


The adaptive grid methodology used in this study is based on the work of Nakahshi and 
Deiwert [ 31 , 32 ]. In this method, the solution adaptive grids are generated based on a tension 
and torsion spring analogy. The tension spring connects adjacent grid points to each other and 
controls the grid spacing so that the grid density is highest in regions where strong gradients 
exist, such as shocks and shear layers. Torsion springs, which are located at each grid node, are 
used to prohibit the grid from becoming excessively skewed. 


B.l One-dimensional Grid Adaption 


Consider a one-dimensional flow in the ^-direction, for which the flows domain has been 
discretized into a grid with constant spacing. Let point A on this grid (Fig. B.l) be connected 
to points B and C by tension springs whose constants are kj t j_i and kj, j respectively .The grid 
points are redistributed along the 77-constant lines by computing the spring constant based on 
the gradient of a given flow parameter, fj, as: 


, _ 1 | C'iI/.j+i ft,j I 

t,j — 1 1 } r~ 

v 5 »,y+i s i,j) 


(B.l) 


where Ci is a constant, and sy is the arc length calculated from point (i,l) along a 77, coordinate 
line. The new location of grid point Sij is computed using the relationship of Eq. (B.l) as: 


ki,j{ 9 i,j + 1 s ij) ki,j—i{ s t,j s i,j- 1) — 0 (B. 2 ) 

If the flow is strickly one-dimensional, the use of tension springs will allow for sufficient control 
over the adaptive grid; however, if the flow is two-dimensional and complex, the grid can 
become excessively skewed (an undesirable effect for most finite volume/difference solvers). To 
alleviate this problem, torsion springs are added to each grid node. The torsion springs control 


63 



the inclination of line DA to that of a reference line. If the torsion spring constant is denoted 
by H ( Fig. B.2), the force at point A is computed as: 


Ftorsion = -Hi-ij($ DA ~ 4>) ( B - 3 ) 

where is an inclination of line DA and <p is the inclination of the reference line. The 
reference line is chosen as an extension of FD to avoid kinks in the £-line at point D. A balance 
equation for the complete spring system is given as: 


s i, j) Hi— <f>i—ij) 0 


(B.4) 


If e and 4> are written in terms of the arc length sy, of the intersection of the reference line DA’ 
and the tj\ coordinate (Fig. B.2) Eq. (B.4) can be rewritten as: 


kij-lSij-l — (ki t j + ki,j - 1 + Hi-l,j) s i,j + ki,j^i,j+l — Hi-ljSij 


(B.5) 


This results in a tridiagonal system of equations for sy which is easily solved. In this analysis, 
only the torsion force on the upstream side ( 77 i_i ) influences the distribution at ?/;. This allows 
for the use of a simple marching type algorithm to be used in the solution of sy. 


B.2 Extension of Grid Adaption to Multi-dimensions 

The basic algorithm described above can be extended to two-direction adaption by performing 
multiple sweeps. In each sweep, the adaption is applied in a single direction only. This is 
analogous to the alternating direction implicit (ADI) logic used in the flow solver described in 
chapter three. This concept is easily shown in Fig. B.3. The two-dimensional adapted grid 
is constructed from two single dimensional adapted grids. Once the grid is “adapted” to the 
flowfield gradients, the solution field is interpolated onto the new grid using second-order, one 
dimensional Lagrange interpolations. 
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Figure B. 1 Tension spring analogy for the one-dimensional adaptive grid 


n i 



Figure B.2 Torsion spring analogy for the adaptive grid 



Figure B.3 Construction of a two-dimensional adapted grid 
using one- dimensional sweeps 
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